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ABSTRACT 



This thesis develops the theory for and presents a 
digital computer program capable of determining the natural 
frequencies of a three-dimensional piping system having 
arbitrary configuration. The analysis uses the method of 
transfer matrices. Piping hangers, loops, and complex 
branches (branches emanating from branches) have not been 
included in the analysis. A distributed mass model is used 
for straight pipe sections, while mass is lumped for curved 
sections. Inclusion of shear deflection and rotary inertia 
is optional. 

Several piping configurations are analyzed using the 
program; the results are compared with analytical solutions 
or values from the literature to demonstrate the accuracy 
and integrity of the program. 
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I. INTRODUCTION 



A. THE PROBLEM 

From an engineering standpoint it is important to be 
able to analyze and predict the vibration characteristics 
of structural elements and assemblies. The mathematical 
aspects of several approaches to such analysis have been 
known for a long time, but the volume of arithmetic opera- 
tions dictated a gross simplification of complicated systems 
in order to reduce the arithmetic to manageable limits. 

The development of the electronic digital computer within 
the last two decades, however, has provided a means for 
the practical and accurate analysis of quite complicated 
vibrating systems without great simplifications which cast 
doubt on the integrity of results. In turn, this has pro- 
vided a stimulus for further developments in theory which 
take advantage of computer capabilities. 

Most methods of analysis of actual systems in which both 
mass and elastic compliance are continuously distributed, 
involve a discretization of one sort or another; usually 
the total mass of the system is considered to be lumped at 
a finite number of distinct points, the masses being con- 
nected by massless spring elements. It is almost universally 
assumed that the elastic characteristics of the system are 
linear and that the excursions of points in the system from 
their equilibrium positions are sufficiently small that no 



6 



geometric nonlinearities are involved. Such an idealization 
leads to n second order ordinary differential equations. 
Provided vibrations are assumed to occur isochronously 
and no damping is present, these equations reduce to a 
system of n algebraic equations in the n amplitudes and 
frequency squared. 

The most convenient and expedient way to cast this 
problem is in matrix form. Thus engineers using these 
methods deal with such quantities as the mass-matrix, 
the stiffness matrix, the modal matrix, etc. 

On the other hand, for certain simple cases of distributed 
mass and elasticity, it is not necessary to lump the mass 
into a finite number of concentrated masses. Instead the 
vibration problem can be formulated as a partial differential 
equation which can be solved directly. Recently, this point 
of view has been extended for the analysis of more complicated 
cases than can be treated by the conventional use of partial 
differential equations. Namely, the system being analyzed 
is divided into a number of sub-systems, a satisfactory 
analysis existing for each one; then, by a procedure which 
will be described in greater detail later, these sub-solutions 
are combined to obtain the solution applicable to the original, 
complicated configuration. This procedure, known as the 
transfer matrix method , evolved from a method described by 
H. Ilolzer fifty years ago for the analysis of torsional 
vibrations. Applications of Holzer's viewpoint to the analysis 
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of flexural systems was made independently by M. A. Prohl 
and N. 0. Myklestad about thirty years ago. The essentials 
of these methods were abstracted by a number of engineers 
and mathematicians who developed and formalized them into 
what is now known as the transfer matrix method. The names 
of S. Falk, K. Marguerre, E. C. Pestel, and F. A. Leckie, 
among others, are associated with this development. A book 
by Pestel and Leckie [Ref. 9] represents the fullest and 
most accessible treatment of the subject available today. 

The transfer matrix method is most naturally applicable 
to the analysis of "chain-like" structures composed of 
elements which are adjoined at distinct points and for 
each of which a satisfactory vibration analysis exists. 

The simplest application is to the case where such a chain 
of elements has but two ends. However, it can be applied 
to structures having slightly greater topological complexity, 
but it simply is not appropriate for the analysis of struc- 
tures which are highly reticulated or in which major elements 
are joined along lines or surfaces rather than at distinct 
points . 

Piping systems are substantially chain-like and topologi- 
cally simple. Accordingly, it was natural that the transfer 
matrix method be applied to the analysis of vibration in 
piping. In this country the first such effort which was 
reported was that done by G. E. Fink in his Naval Postgraduate 
School thesis of 1964. At about the same time similar 
efforts were being made in Japan under the sponsorship of 
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the Japan Society of Mechanical Engineers. It is known 
that the Japanese have produced useful digital computer 
programs for the analysis of piping using the transfer 
matrix method, but it has not been possible to learn any 
of the details since these programs remain under the pro- 
prietary control of the Japanese government. 

Fink developed a program which has convenient generality 
and gives good results, but which is limited to systems 
lying entirely in a single plane. His program is capable of 
making separate determinations of "in-plane" and "out-of- 
plane" vibrations, which, for the case of uniplanar con- 
figurations, are uncoupled. Further work on the theory 
was done by W. S. Baird, Jr. and J. L. Simmons in their 
Naval Postgraduate School theses; this work mainly related 
to systems having appreciable topological complexity. In 
his Naval Postgraduate School thesis, Y. S. Kim returned 
to the development of Fink's program, adding several useful 
features, including the incorporation of several alternative 
mathematical techniques for accurate solution. Some of these 
techniques reduced solution time below that required in Fink's 
original program; however, Kim's program is still limited to 
uniplanar systems. 

Accordingly, the writer undertook to develop a program 
capable of treating truly three-dimensional configurations. 

The remainder of this thesis presents the analytical back- 
ground of this development, lists the program which includes 
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rules for its operation, and gives evidence concerning the 
high accuracy of the results. 

It should be made clear that no claims are made concerning 
the relative usefulness of the transfer matrix method as 
compared to lumped-mass methods which are presently in 
widespread use for analysis of piping vibrations. In 
general, the writer believes that the results which are 
obtained by the transfer matrix procedure are substantially 
more accurate than can be obtained in most cases using lumped- 
mass analysis. However, the delicacy of the calculations 
involved in the transfer matrix method seems to imply 
that only a limited number of the lowest frequencies (and 
their associated modes) can be obtained using present pro- 
gramming methods and computer hardware. Greater word- length 
in the computer, or availability of fast multiple precision 
capability, would permit obtaining a larger number of 
frequencies and modes via the transfer matrix method. 

Briefly, the motivation for the development reported 
herein has been simply that the transfer matrix method is 
clearly appropriate for the analysis of piping vibrations 
and that the application should be made. It seems likely 
that in some cases and for some purposes the transfer matrix 
method may be more convenient and economical while in other 
cases lumped-mass methods may be preferable. 
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B. THE TRANSFER MATRIX METHOD 



Specific details of the application of transfer matrix 
methodology to the analysis of piping vibrations will be 
given in subsequent parts of this thesis. Background 
material and a general exposition of the method is most 
readily available in the book by Pestel and Leckie [Ref. 9] 
and it would be pointless to include herein a paraphrase 
of such material. However, for the reader who is not familiar 
with the ideas behind the transfer matrix method, a very 
brief exposition is given in this section. 

An appropriate "state vector" is defined which presents 
force-type and deflection- type information capable of specify- 
int the "state" at each of several points in the configura- 
tion at which sub-elements are joined together. These 
quantities specify the configuration and internal force 
system in the structure at its extreme isochronous, deflected 
position. For example, in a simple torsional system, there 
are only two quantities in the state vector, the angle of 
rotation from the equilibrium position and the torque 
transmitted from one element to the next. The elastic 
compliance and the inertial characteristics of an element 
permit obtaining equations relating the quantities in the 
state vector at the left end of the element to those in 
the state vector at the right end of the element. When 
these relations are put in matrix form it is found that 
the state vector at the right can be written as the product 
of the state vector at the left premultiplied by a square 
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matrix, called a transfer matrix, the elements of which 
represent the elastic and inertial properties of the element. 

In applying the method, a state vector is chosen for the 
left end of the assembly of elements in such a way as to 
satisfy conditions of restraint at the left. Then, by suc- 
cessive multiplication by transfer matrices representing 
individual elements encountered in order, progressing from 
the left end of the assembly toward the right, one arrives 
at a representation for the state vector at the right end. 

The individual transfer matrices are functions of isochronous 
frequency; accordingly the right end matrix is a function of 
frequency. When frequency is adjusted, any particular 
frequency which permits satisfaction of restraint at the 
right end is a natural frequency of the system. The several 
state vectors at the junction points give the deflected 
configuration, or mode shape, corresponding to this frequency, 
when this value is substituted in the expressions for the 
state vectors. 

For simple cases it is possible to carry out the multi- 
Pl ications in literal form, keeping the frequency as an 
unknown parameter. Then the right end restraint conditions 
provide a polynomial equation in the frequency which can 
be solved by customary methods. However, in relatively 
complicated cases it is not feasible to retain the frequency 
as a literal parameter. In these cases a definite numerical 
value is assumed for frequency and this leads to a numerical 
measure of the failure to satisfy right end restraint 
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conditions. By assuming different values of the frequency 
and by constructing a curve of measure-of-f ailure versus 
frequency, the correct frequencies can be obtained. 

The transfer matrices which describe the individual 
elements may be a theoretically complete and full repre- 
sentation, or they may themselves by approximate representa- 
tions. In the application to piping systems it is convenient 
to represent a straight length of pipe by a transfer matrix 
which is theoretically complete and exact. For systems 
composed of only straight lengths of pipe there is no 
approximation whatsoever in the transfer matrix method 
(other than that implied by internal round-off in the com- 
puter) ; this may be compared with the fact that all lumped- 
mass procedures definitely imply an approximation and thus 
result in some theoretical error. 

Curved pipe elements (elbows and bends) could be treated 
exactly since adequate theory exists for the construction 
of theoretically exact transfer matrices with which to 
represent them. However, there are severe computational 
difficulties involved with such exact representation. 
Accordingly, with the awareness that the total length of 
curved pipe is but a small fraction of the total length 
of both straight and curved pipe in the average configura- 
tion, it was decided to represent curved pipe by a succes- 
sion of massless curved sections having proper elastic 
properties, alternating with appropriately located point 
masses. For each of these idealized elements, it is 



possible to obtain the corresponding transfer matrices with 
no computational difficulties. Thus, elbows and bends are 
represented in much the same way that is done with the 
lumped-mass methods mentioned previously. It is believed 
that the approximation introduced by this point-mass 
representation of a small part of a configuration is quite 
tolerable for practical engineering purposes. However, it 
is certainly possible to modify the program reported herein 
so as to provide for exact representation for curved ele- 
ments. This would not involve major changes but only the 
development of a subroutine to replace those given herein 
for curved elements. 
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II. PROGRAM DEVELOPMENT 



A. THE STATE VECTOR, TRANSFER MATRIX, AND FREQUENCY 
DETERMINANT 

1 . State Vector 

The state vector at point i of an elastic system 
is a column vector whose components are the generalized 
displacements, D^ , and the generalized forces, F^, at that 
point. In matrix notation, we have 




For a beam free to move in three dimensions with 



longitudinal, torsional and flexural elasticity, the state 



vector takes the form 
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deflection in x-direction 
deflection in y-direction 
deflection in z-direction 
rotation about x - axis 
rotation about y-axis 
rotation about z-axis 
axial force in x-direction 
shear force in y-direction 
shear force in z-direction 
Torque about x-axis 
Moment about y-axis 
Moment about z-axis 



The order in which the individual displacements and forces 



occur in the state vector is of no particular importance 



so long as consistency is maintained. Program VIBREL 
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(for VIBRation of ELastic piping) , the program developed 
herein, employs the arrangement given above. 

2 . Sign Convention and Coordinate System 

A right-handed cartesian coordinate system is used 
to describe the local generalized forces and displacements, 
the x-axis coinciding with the centroidal axis of the pipe 
as shown in Fig. 2.1. 




Figure 2.1: Vectorial Representation of State Vector 

Displacement and Forces 

A cut made at any location of the pipe will expose two 
faces; the face with its outward normal in the positive x-axis 
direction is considered positive. The state vector parameters 
are positive if their vector representations coincide with 
the directions prescribed by the cartesian coordinate 
system at the positive face as in Fig. 2-1. Rotation and 
moment vectors are portrayed according to the right-hand- 
screw rule by a double line with an arrowhead. 

3 . Transfer Matrices 

The transfer matrix is defined as the matrix which 
relates two state vectors at successive selected points of 
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division in the elastic system. 

z . = U . z . . 

1 11-I 

where Ih is the transfer matrix. The transfer matrix must 
be square, and in the case of a spatial system must have 12 
rows and 12 columns. When the transfer matrix relates the 
flexibility between the points i and i-1 of a continuous 
system it will hereafter be referred to as a field matrix, V. 
When a discontinuity arises in either force or displacement, 
as with a point mass in a lumped system, the transfer matrix 
relating the state vectors on either side of the discontinu- 
ity is called a point matrix, P. The term "transfer matrix" 
will henceforth be used with regard to the matrix, U, which 
includes mass and flexibility in the relationship between 
the state vectors at two distinct points. A complete discus- 
sion of transfer matrices, with examples and derivations, can 
be found in Ref. 9. 

4 . The Elimination Process and Frequency Determinant 

The process of combining the system transfer matrices 
to find natural frequencies can be demonstrated by a simple 
planar example using the Myklestad method of replacing a 
distributed mass beam with lumped masses connected by flexi- 
ble elements (Fig. 2.2). Assuming that the field and point 




i 



*m 



£ <7 



i 



Figure 2.2: Lumped Mass Representation of 

a Planar Beam 
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matrices for a massless flexible beam and a point mass are 
already known, then 



L R n L j -.r R 
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where z . = A 
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. and superscripts L and R indicate left side 



or right side, respectively. Combining the above yields the 
relationship 



z 0 = V~PV,z = Uz 
2 2 1 o o 



( 2 - 1 ) 



where U is the system transfer matrix. Written in full 
this relationship is 
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the a^. being vibration functions of mass, stiffness, and 
frequency. Application of the boundary conditions w q = = 

\p o = \p 2 = 0 reduces the matrix equation to the following 
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Because the elements of this square matrix are functions of 
frequency, its determinant is called the frequency determi- 
nant of the system. Note that this 2X2 matrix is a 
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submatrix of the 4X4 matrix appearing above. A non- 
trivial solution of the equations 2-2 requires that the value 
of the frequency determinant (A) equal zero. A value of the 
frequency, co, which causes A to vanish is called a natural 
frequency of the system. 

The natural frequencies could be solved for directly 
provided the expressions for the a^ were kept in terms of 
the circular frequency, w, but in practice these expressions 
are forms much too unwieldy to be handled explicitly. 

Program VIBREL evaluates the frequency determinant of 
order 6X6 for successive numerical values of frequency, 
zero values of the frequency determinant corresponding to 
natural frequencies of the system. 

Let the state vector be composed of 2r quantities. Then 
to compute U in Eq. 2-1, two multiplications of 2r X 2r 
matrices are required. Noticing that application of the 
boundary conditions to the state vector at point 0 will 
cause r of its elements to be zero, we find that r columns 
of U are multiplied by zero, and thus play no part in the 
calculation. Kim [Ref. 6] introduced the concept of the 
"state matrix,” S, such that the left end state vector, z q , 
of 2r quantities can be written in terms of a compressed 
state vector, z q * , having only r quantities, viz., 

z = Sz * 
o o 

where S is the state matrix at the left end, having r 
columns and 2r rows, its elements consisting only of zeros 
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and ones. All elements of the state matrix are zero except 
that in each column, unity appears in the same row position 
as the corresponding non-zero term in the state vector, 
boundary conditions having been applied. If the state 
matrix is used in the multiplication scheme, then we can 
multiply 2r X r matrices by 2r X 2r matrices and reduce 
the calculation considerably. The resulting 2r X r system 
transfer matrix requires only the application of right end 
boundary conditions to yield the frequency condition. 

To illustrate this state matrix concept most simply, 
consider again a two-dimensional case in which the state 
vector consists of four quantities. The three dimensional 
case is an obvious generalization of the following discus- 
sion. Taking strictly the state vector approach we can 




Figure 2.3: Single Branch Planar System 



express the relation between the state vector at 2 and 
that at 0 by 



z 7 = U,P,U.,z rt 
2 2 1 1 o 



(2-3) 



where and are straight section transfer matrices and 
the point matrix P accounts for the effects of the branch. 
Carrying out the matrix multiplication and exhibiting the 
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system transfer matrix and state vectors which consider 
bending only, we have 
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The boundary conditions have been applied to the left end 
state vector; the quantities and A 2 represent the unknown 
quantities of and respectively. Application of the 
right end boundary conditions yields a frequency condition 
of the same form as Eq. 2-2. 

Using the state matrix approach instead, Eq. 2-3 
becomes 



z 9 = U-P.U^Sz * 
2 2 1 1 o 



(2-5) 



where S is the state matrix 



z * contains the non- 
o 



0 0 
0 0 
1 0 
0 1 

zero quantities of the left end state vector. After multi- 
plication of the matrices of Eq. 2-5, the number of individual 
multiplications being considerably less than before, the 
following relation is obtained: 



ft* - 

w 




* a 

13 


...i 


< 


^ = 


a 

23 


a 24 


y 




a 


a 


Y z 




33 


34 


' M 

z 


0 


a 43 


a 44 


L J 


z 


- 


- 



r *■* 

A i 

A. 



(2-6) 



21 



We see from this result that inclusion of the state 



matrix compresses the left end state vector to r quantities 
while providing the same useful information as the matrix 
of Eq. 2-4, inasmuch as the same frequency condition is 
obtained after the right end boundary conditions are 
applied. Also, if the relative magnitudes of the unknowns 
and have been determined, the mode deflections and 
forces at a point can be calculated simply by multiplica- 
tion of the compressed left end state vector, consisting 
of r quantities, by the 2r X r system transfer matrix up 
to that point. The appropriate mode frequency is used in 
the computation. This will be discussed in greater detail 
in Chapter 3. 

B. 3-D PIPING GEOMETRY 

1 . Working Points and Piping Nomenclature 

The geometry of the piping system as it is used in 
program VIBREL is based on working points. A working point 
of the system occurs where cross section or properties 
change, where branches join the main member, where the 
piping changes direction, and at the ends of the main member 
and branches. These points are measured in cartesian coordi- 
nates from an origin located arbitrarily but fixed in space. 
For purposes of this discussion we shall distinguish between 
the two types of working points at which piping direction 
changes. At a corner point (Fig. 2.4) the piping direction 
changes with no discernible radius of curvature. At a 
bend point (Fig. 2.5), which is defined as the intersection 



of the tangents to the two ends of a curved pipe, the radius 
of curvature is finite. 




Figure 2.4: Piping Corner 




The main member is designated as the longest continuous 
run of pipe from which all branches emanate. Program 
VIBREL cannot handle branches emanating from branches. 

2 . Unit Vectors and Coordinate Transformations 

If a set of mutually orthogonal x, y, and z unit 
vectors are formed at the main member and branch starting 
points, and prior to and following every every directional 
change, we can determine, by angular differences in succes- 
sive sets of unit vectors, the appropriate coordinate trans- 
formation angles. The local coordinate systems described 
by these unit vectors are formed with the piping system 
in its quiescent configuration. 

Let us take, for example, the two sections of pipe 
pictured in Fig. 2.6 determined by the three working points 
E, F, G, where F is a corner point. Also let the two sections 
of pipe, EF and FG, be represented by vectors X and S, the 
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Figure 2.6: Vector Representation of Directional Change 



components of which have been determined from the coordi- 
nates of E, F, and G. Then the angle y can be determined by 
the following relation: 



r~ 




The sign will be determined later. Since the vectorial 
representation of state vector quantities in Fig. 2.1 
requires that the x-axis coincide with the centroidal 
axis of the pipe, let 



X I A | 

For the purposes of succeding derivations we desire that 
e^ be in the plane of A and B and point outward from the 
directional change. Hence 



B x A 
| B x A | 



and 



e = e x e 
y z x 



By the right hand screw convention the rotation of amount 
Y will always be in the negative e y sense. The set of unit 
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vectors, U^, U , U z are found simply by rotating the initial 

unit vectors about the local z-axis by the angle y. 

Now consider the difference in orientation between 

the 0 , U„, 0, vectors and those prior to the next directional 
x y z 1 

change farther down the pipe,U x ,, 0 , , and U z , . These have been 

constructed in a manner similar to e , e , and e . Referring 

x y z ° 

to Fig. 2.7, the angular difference between the 




Figure 2.7: Unit Vector Rotation about x-axis 

two adjacent sets of unit vectors is merely the rotation a 
about the local x-axis. Alpha may be determined by the 
expression 

M - COS (Oy-Uy'J 

This rotation is in the positive (negative) 0 sense accord- 

A 

ing to whether the triple product (U^xU ,) *0^’ is positive 
(negative) . 

In this way, the state vector is always expressed 
in the appropriate local coordinate system. For the simple 
system of Fig. 2.8 




Figure 2.8: Planar System with Directional Change 



z, = U, z 
1 1 o 
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IT L 

U z, 
Y 1 



R 



z 2 = U 2 Z 1 



which combining gives 



z, = U,U U,z 
2 2 y 1 o 

where is the transformation matrix corresponding 
rotation y about the z-axis. For a 3-D system with 
vector defined as in section 2.1 this matrix is 



to a 

the state 



U = DIAG [T , T , T , T ] 
Y *- y 7 y "Y* Y J 
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-SIN y COS y 0 

0 0 1 



Similarly, for the x-axis rotations 



U = DIAG [T , T , 
a L a ’ a ’ 





where T 
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10 0 
0 COS a SIN a 
0 -SIN a COS a 



C. 3-D TRANSFER MATRICES 

1 . The Straight Pipe Transfer Matrix 

In order to maximize accuracy, a distributed mass 
model was decided upon for straight pieces of pipe. The 
program uses the straight bar transfer matrix cited in the 
appendix of Ref. 9, but revised to conform to the ordering 
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of the state vector described herein. The revised matrix, 
which appears in Appendix B, combines the vibration func- 
tions for longitudinal, torsional, and flexural modes and 
can include or neglect rotary inertia and shear deflection. 

2 . Development of the Curved Pipe Field Matrix 

Due to the complexity and increased machine time 
involved in consideration of a distributed mass model for 
curved sections of pipe, a lumped mass treatment was chosen. 
This, coupled with the distributed mass approach for straight 
pipe, yields staisfactory accuracy combined with minimal 
computer time. 

Given the compliance of a piece of massless curved 

pipe at its center of curvature, we seek the field matrix 

D. 

V such that z^ = Vz^ (see Fig. 2.9) where z^ = { p 1 } . The 




Figure 2.9: Massless Curved Pipe 

compliance matrix is defined by the relationship 
D = C F where D is a vector of generalized displacements 
observed at P and F^ is some vector of generalized forces 
applied at- P. Point P is assumed to be attached by a rigid' 
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bar to one end of the pipe while the other end of the pipe 
remains fixed. The compliance matrix for a pipe bend and 
its derivation are presented in a paper by J. E. Brock 

[Ref. 2]. 

The analysis requires that we find the compliance 
matrix at point L which is obtained by the appropriate 
rotation and translation of coordinates from point P such 
that 

C L = B PL L PL C P L PL B PL 
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for rotation. 
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The matrix equation for the generalized displacements at 
point L can now be written. 



d l c l f l + b rl l lr d r 

where the rotation and translation matrices use the above 
notation. Solving for D p we arrive at 
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The equation for the generalized forces at the right end 
of the pipe has the form 



f r " l rl b rl f l 



(2-8) 



If equations 2-7 and 2-8 are combined in matrix form, we 
have the result 
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By the sign convention discussed in section 2.1, point L 
is located at a negative face, so after appropriate sign 
changes are taken into account, the field matrix for the 
bend becomes 

T Tp T j T T n T r 

LR LR ' LR LR L 

v = r 

o ! " l rl b rl 

w i 

This matrix in its explicit form, as used in program 
VIBREL, appears in Appendix B. 

3 . Development of Curved Pipe Point Matrix 

The mass of a pipe bend or elbow is taken into 
account through the transfer matrix procedure by lumping it 
at the center of mass and including its effect on the system 
at point 0 (Fig. 2.10), the left end of the elbow section. 
The resulting point matrix relates only the forces across 
point 0 since the deflections are continuous. 
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Figure 2-10: Curved Pipe Showing Location 

of Center of Mass 
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A force summation at point 0 yields the following 
scalar equations (the quantities appearing in this develop- 
ment are defined in Appendix A, subscripts L and R refer 
to the left and right of point 0, respectively), 

N r = n l - mw 2 (u - y0) (2-9) 

V yR = V yL " mw2 ( v + * 6 ) (2-10) 

V zR = V zL ~ m(Jj2 ( w ' Y$ " (2-11) 



Similarly, a moment summation at the same point produces 
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Combining equations 2-9 through 2-14 
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and the resulting point matrix equation is: 
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An inertia tensor analysis was used in the determination 
of the mass moments of inertia about point 0. It can 
readily be shown that the inertia tensor at point 0 is 



,-2 



?o - ^r-t (J+SIN2Y) n _ sin 2 y (iljii) + (Y. SIN2Y) j ■ j 
+2pAR 3 [y-SINy]u + pAR 3 [(-£ + SIN |^)ii 
+ (1 - COSy - + - (p- - 2SINy + j j ] 



where it is assumed that pipe wall thickness is small 
compared to pipe radius. (The symbols used here are defined 

— -~y 

in Appendix A.) Then I = i. I • l, I = j • I • j, 

J xo o ’ yo J o J ’ 

— = — -y -y 

and I = k • I • k, where l, i , and k are unit vectors 

along the local coordinate axes at point 0. The ensuing 

moments of inertia and point matrix are presented in their 
entirety in Appendix B. 



D. BRANCHED SYSTEMS 

Any branch joining a piping system has a significant 
effect on its dynamic behavior. While the deflections 
across the junction point are continuous, there is a dis- 
continuity in the forces, the magnitude of which is 
dependent on the displacement at the junction point and 
the nature of the branch. 
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Consider the simple branch system of Fig. 2.11. 




Figure 2.11: Unit Vector Orientation at a Branch Junction 



Using the state matrix concept discussed in section 2.1 we 
have the following relationship between the state vector at 
point D and the state vector at point B (the bar denotes 
branch coordinate scheme) . 
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The zero quantities of the state vector at D have been 
eliminated by application of the boundary conditions and 
introduction of the state matrix, and the branch transfer 
matrix has been partitioned into two submatrices and 
The local coordinate system of the branch at point B is 
represented by the unprimed set of unit vectors. From the 
relationship 2-15 it can be seen that Dg = R-^Zg and 
Fg = R 2 z p> so that when Zg is eliminated, the equation 
between the forces and displacements at point B in the 
unprimed coordinate system becomes Fg = P^R^ (2-16). 



The vectors of the generalized forces and displacements 
in the branch system must now be transferred to the main 
member system represented by the unit vectors s in Fig. 2.11. 
This can be achieved by one z-axis and two x-axis transfor- 
mations of the unprimed coordinates. We first form the 
primed set of unit vectors in a manner identical to that for 
a corner point (section 2.2), the three working points 
describing the corner being E, B and C. An x-axis rotation 
a^ will align the unprimed with the primed coordinates. 

t 

Next, a rotation y about the z -axis will allow u , to 
coincide with the centoidal axis of section BC resulting in 
the double-primed coordinate scheme. Finally, rotation 
through the angle a about the x"-axis will align the branch 
system with the main member coordinates, s. The total trans- 
formation of the displacements can be expressed by 
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similarly for G , and 
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Likewise for the forces, Fg = GFg (2-18). Rearranging 



relations 2-17 and 2-18 and substituting in 2-16 yields 
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* 




Noting that G 
equation 



G'\ = «2 R ; 1g ' 1d B • 



- 1 T 

= G we now form the force displacement 
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Referring to the free body diagram in Fig. 2.12 we find 
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Figure 2.12: Generalized Forces at a Branch Junction 

that the main member force equation at point B can be written 



F R = F L + F 



B 



I - 1 T T T 

= F + GR-R, G ,G G D„ 

21 al y a B . 



,R 



Thus, since = D = D , the point matrix relating the state 
vectors to the left and right of a branch junction point of 
a 3-D system is given by 
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The foregoing development assumes that there is no 
curvature of the main member or branch at point B and that 
no abrupt change of direction of the main member occurs at 
that immediate location. 
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For two or more branches joining the main member at 
a single point, the system transfer matrix just prior to 
the branch point is simply multiplied in turn by the point 
matrix for each of the branches. 

E. COMPUTER IMPLEMENTATION OF TRANSFER MATRICES 

It would be an extremely tedious job to attempt explicit 
assemblage of the elements of the system transfer matrix in 
terms of the circular frequency. Subsequent solution for the 
roots of the expanded frequency determinant would also be 
impractical for all but the simplest of systems. It is 
for these reasons that digital computation becomes a 
necessity . 

Program VIBREL, which incorporates the theory of 
preceding sections, was developed to utilize the speed with 
which the complete system transfer matrix can be formed time 
after time for various values of frequency. A typical graph 
of frequency determinant A (to) versus circular frequency to 
is shown in Fig. 2.13. The values of to at which the graph 




Figure 2.13: Graph of Frequency Determinant vs. Frequency 
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crosses the horizontal axis, namely those values of w for 
which A(w) = 0, are the natural frequencies of the system. 

The program, presented in Appendices C through E, was 
coded using Fortran language [Ref. 7] in double precision 
to minimize the effects of round-off and truncation errors. 
Since the prime objective was to produce a working program 
of demonstrable accuracy and reliability, limited attention 
was focused on the niceties of programming intended to 
reduce execution time. However, Kim's state matrix concept 
is included and this effects a substantial saving of time 
as compared to other procedures. There is no reason to 
believe that the actual programming is particularly ineffi- 
cient; however, it is not unlikely that close scrutiny of 
program details might point oirt some places where slight 
revisions would effect economies. 

The transfer matrices used in computation of system 
frequencies are those presented in Appendix B. 

Input formats were designed to provide satisfactory 
simplicity for the user in preparing data decks, while 
output formats were chosen to check for correct input as well 
as to provide a detailed summary of results. The I/O 
formats are discussed more thoroughly in Appendix E. 
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III. DISCUSSION 



A. SCOPE OF THE SOLUTION 

In its present form, program VIBREL is capable of 
analyzing a piping system having few restrictions regarding 
topology. The system must, however, be composed of a 
main member from which can emanate as many as fifty branches. 
The program as listed in Appendix E can accommodate no more 
than two branches joining the main member at a particular 
point. In addition, there can be no branches emanating 
from other branches, nor can there by any curvature of the 
main member or branch at a branch junction point. 

The total number of sections capable of being analyzed, 
presently set at one hundred, can be increased by follow- 
ing the procedure described in Appendix C. Likewise the 
number of branches can be augmented from its present limit 
of fifty. For each increase of one section, eighty addi- 
tional bytes of computer storage are needed above the 140,000 
bytes now required for the entire program. For most computer 
installations this allows considerable leeway in the size 
of the piping system which can be handled by the program. 

To provide an idea of the computer time involved in 
the frequency analysis, the twelve- section piping system 
of the sample problem in Appendix G required about four 
minutes for computation of four modes and almost fifteen 
minutes for eighteen modes. The IBM 360 computer compiles 
the program in thirty-nine seconds. 
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The usual assumptions of linearity as well as homogeneity 
and isotropy of the structural material have been incorporated 
in the analysis. All sections of the system are assumed to 
vibrate isochronously with negligible damping. 

B. NATURE OF ERRORS AND INACCURACIES 
1 . Round-Off Errors 

As piping system complexity increases the number 

of multiplications required to form the system transfer 

matrix grows proportionately. Each time a multiplication 

occurs, there is some round-off and loss of accuracy; this 

is a function of the significant figure capacity of the computer. 

FORTRAN, in conjunction with the IBM 360 and double 

precision arithmetic, has the capability of carrying numbers 

-78 

of 16 significant digits with an exponent range of 10 to 
7 8 

10 . Simpler systems can be checked for accuracy and it is 

evident that this significant digit capacity keeps round- 
off errors negligible. Since with more complex systems 
other forms of errors creep into the results, we have no 
way of predicting the exact contribution of round-off 
errors to total inaccuracy of frequency. 

It is reasonable to assume that round-off plays 
a more significant role in the 3-D piping system accuracy than 
in a planar case with equivalent sections because of the 
greater number of coordinate transformations that are neces- 
sary for system description. 
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2 . Lumped Mass Idealization 



For sections of curved pipe, lumped mass connected 
by massless springs is used as a model for the real system. 
This idealization leads to inherent inaccuracy in the analy- 
sis of systems containing bends. 

3 . Zeros of the Frequency Determinant 

Once it has been established that the frequency 
determinant has crossed the zero axis, Program VIBREL uses 
Mueller’s method of successive bisection and inverse parabolic 
interpolation to determine the zero. Care should be exer- 
cised in considering as significant only the same number of 
places as in the acceptability criterion which was specified 
for the solution. (See Appendix C.) 

4 . Other Sources of Inaccuracies 

The imperfect nature of the real system and inevitable 
errors in the measurement of working points will render the 
computer solution only a good engineering approximation of 
What can be expected in actuality. In the discussion of 
accuracy in succeeding sections, inaccuracies caused by real 
system imperfections will not be considered. 

C. ACCURACY AND INTEGRITY OF SOLUTION 

1 . Solution Accuracy 

Confidence in the accuracy of the program solution 
was established by three comparison tests, the results of 
which were tabulated in Appendix F. 
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a. Analytical Comparison 

For single straight sections of pipe a closed 
form solution was obtained from the governing fourth order 
differential equation with the appropriate fixed, free, or 
ball joint end conditions. The largest variation between 
the natural frequencies calculated analytically and those 
computed by program VIBREL was 0.009%. 

b. Reference Value Comparison 

Where the exact closed form solution was not 
available, recourse to the literature provided some compari- 
son frequencies. For the simple curved pipe shown in 
Appendix F, Refs. 5, 6, and 8 contained a range of frequency 
values for the first two modes. VIBREL results for the 
same curved pipe fell within this range in each of the two 
modes examined, with an average variation from the reference 
values of about 4%. Because, for a curved pipe, program 
VIBREL includes bend flexibility factor and an accurate 
value for shear distribution factor taken from Ref 2, the 
values computed are considered to be at least as accurate as 
those contained in the literature. 

For a single branch system composed of straight 
pipes only, comparison frequencies for three modes were 
available from Ref. 6. The largest difference between 
reference and VIBREL in this case was 0.06%. 

c. Dual Analysis Comparison 

Dual analysis is based on the fact that the system 
transfer matrix is directionally dependent. One end of the 
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system main member is designated as the starting point; then 
the component transfer matrices are multiplied together in 
sequence proceeding toward the other end. If an error 
were to occur in the structure of one of the transfer 
matrices or in the method of their combination to form the 
system transfer matrix, the results from starting at opposing 
ends could differ significantly. 

An improperly constructed state matrix would 
also be evident from a dual analysis. 

The straight pipe systems and branched systems 
showed no significant differences when subjected to dual 
analysis. The maximum difference observed in eighteen 
modes of the complex system dual analysis was 2.67%, 
while the average difference was 0.58%. 

2 . Solution Integrity 

Solution integrity was established by observing that 
the mode frequencies detected by the program were, in fact, 
those of the system with no omissions. This was accomplished 
by comparison with analytical and reference results. Dual 
analysis comparison was also used in determining that the 
same modes had been detected in either direction. No dis- 
crepancies were noted in any of the systems checked on the 
basis of these comparisons. 

Program VIBREL includes in its output a graph of 
frequency determinant versus frequency in order to check 
whether the frequency increment is such that a mode or 
modes have been skipped in the scanning process. In the 
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first analysis of a branched system which appears in 
Appendix F, the trend of the output plot showed that two 
natural frequencies in close proximity had been overlooked. 
When the program was rerun with the starting frequency incre- 
ment decreased by a factor of 2, these two frequencies were 
detected . 

The maximum number of natural frequencies which can 
be found using VIBREL is related to the significant figure 
capacity of the computer and the physical configuration 
of the piping system. As pointed out by Kim [Ref. 6] , 
at the higher frequencies the columns of the frequency 
determinant approach the point of being parallel; hence the 
numerical value of the determinant includes differences of 
large numbers. Systems with few components generally 
exhibit this parallelism at a lower frequency than the more 
complex ones; however, there is no way of determining before- 
hand where numerical difficulties will be encountered. When 
parallelism is approached, a scattering of the values of 
frequency determinant will be noticed in the output plot. 

Although the exponent limit on the IBM 360 computer 
used for testing the program is ^v78, program VIBREL has been 
coded to cease computation for a particular problem and print 
a message when the exponent of frequency determinant exceeds 
+ 60. 
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D. CONCLUSIONS AND RECOMMENDATIONS 



Based on the results of the accuracy and integrity 
analyses previously discussed and through comparison with 
similar programs developed for planar cases, namely those 
of Fink [Ref 5] and Kim [Ref. 6] , program VIBREL could be 
employed in its present form as a working tool for engineers 
in the dynamic analysis of spatial piping systems. The 
fact that eighteen mode frequencies were obtained for the 
complex system exhibited in Appendix F implies that at 
least the first few and usually the most useful frequencies 
can be determined for practical piping configuration with 
an accuracy consistent with engineering design. 

Minor modifications for adapting the program to a 
particular situation can be made using the guidelines of 
Appendix C. With regard to major revisions, program VIBREL 
was developed with solution accuracy and integrity of first 
importance. It is recommended for future users that changes 
in the coding be made for a minimum computer execution time 
if that becomes a prime requisite. Avenues of investigation 
regarding numerical difficulties at higher frequencies are 
also open to the future user. A similar analysis for the 
planar case was undertaken by Kim [Ref. 6], who utilized 
alternatives to the frequency determinant approach. 

Although the calculation of mode shapes of the piping 
system has not been included in the program because of time 
limitations, an outline of how to incorporate this feature 
is discussed in the following paragraphs. This modification 
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to the program is recommended as a further important step 
in the dynamic analysis of 3-D piping systems. 

Let us consider again, for simplicity, the two-dimensional 
system of Fig. 2.3 which is repeated below. We found that 
the total system frequency condition could be expressed as: 
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(3-1) 




Figure 2.3: Single Branch Planar System 

A vanishing determinant of the coefficient matrix yields a 

natural frequency. At this frequency there is a non- trivial 

solution to Eqs. 3-1 and A^/A 2 can be found. Normalizing 

A^ and gives the modal values of the non-zero quantities 

of the left end state vector, A and A^* . A suitable nor- 

2 2 

malizing condition would be A^ + A 2 = 1. With these values 
known we proceed along the main member calculating the modal 
deflections and forces at each divisional point. This is 
done simply by multiplying the normalized and compressed 
left end state vector by the system transfer matrix (2r x r) . 
In Fig. 2.3 the modal forces and deflections to the left of 
point 1 would be 
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In order to permit unequivocal interpretation of program 
output, it would be advisable to incorporate a point number- 
ing scheme for each point at which modal forces and deflec- 
tions are desired. In addition, these forces and deflections 
will be expressed in the local coordinate systems and a method 
should be incorporated for transforming them to the global 
sys tern . 

When a branch point, such as 1 in Fig. 2.3, is encountered, 
the branch transfer matrix up to point 1 is formed. Thus 
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where and V>^ are the unknown quantities of the state 
vector at the remote end of the branch. Since we know that 
the deflections are continuous across the junction point 



w 






< 


p- = " 




-h 


.B 


>i 



(3-4) 



46 



The quantities at 1L are given in terms of the known values 
A-^* and A£* so that the values of B^ and B ^ which are intrin- 
sic to the IB deflections may be solved for through the 
relations 3-4. The known values of B^ and ma Y now be 
used to find modal forces and deflections in a manner simi- 
lar to that used for the main member. 

Just to the right of the junction point, 1R, the modal 
deflections are the same as at the left; however, the forces 
must be calculated by adding the 1L forces to the IB forces 
in the global coordinate system. 



47 



APPENDIX A 



NOTATION AND NOMENCLATURE 



The following matrix notation is used in the text: 

[ ] - matrix 

{ } - n x 1 column vector 

DIAG { } - square matrix having the given elements 
in its principal diagonal and zero 

elsewhere 

T 

- superscript denoting the transpose of a given 
matrix 

The following vector notation is used in the text: 

A - vector A 
| A | - length of vector A 
A-B - dot product of vectors A and B 
A x B - cross product of vectors A and B 
Table A.l first lists the symbols used in either the 
test or Appendix B, with the corresponding computer program 
variable identifiers. Following these appear the program 
variable names which have no counterpart in the text. 



TABLE A.l 



TEXT PROGRAM 
SYMBOL VARIABLE 



DESCRIPTION 



A 



AREA 



cross sectional area of pipe 

first of two vectors defining a piping 
directional change 



A 



A 
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TABLE A.l (Continued) 



TEXT PROGRAM 
SYMBOL VARIABLE 

B 

B B 



C 

D 

®x UX 

E EY 

F 

G G 

G G2 

a 

G G3 

Y 

G , G1 

al 

i 

I x DIX 



! o 

V 

i 

5 

J 

k 

k 

l 

L 

m 



DJ 



BFF 

DL 



DM 



DESCRIPTION 
translation matrix 

second of two vectors defining a piping 
directional change 

compliance matrix 

vector of generalized displacements 
unit vector in local x-direction 
Young’s modulus 
vector of generalized forces 
shear modulus 

branch coordinate transformation matrix 

branch coordinate transformation matrix 

branch coordinate transformation matrix 

unit vector in local x-direction 

mass moment of inertia of a curved pipe 
about x-axis 

inei;tia tensor about point 0 
6x6 identity matrix 

radius of gyration of pipe cross section 

unit vector in local y-direction 

moment of inertia of pipe section 

unit vector in local z-direction 

bend flexibility factor 

length of a straight section of pipe 

rotation matrix 

mass of a section of pipe 
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TABLE A.l (Continued) 



TEXT 

SYMBOL 


PROGRAM 

VARIABLE 


DESCRIPTION 


M 

y 




moment about local y-axis 


M 

z 




moment about local z-axis 


N 




axial force in local x-direction 


°6 




6x6 null matrix 


P 




point matrix 


R 


R 


radius of curvature of pipe bend 


r 


RBAR 


mean radius of pipe section 


R 1 


R1 


submatrix of branch transfer matrix 


R 2 


R2 


submatrix of branch transfer matrix 


r G 


RG 


radius of center of mass of a curved 
section of pipe 


s 

X 


SX 


main member unit vector just prior t 
branch junction 


t 


T 


pipe wall thickness 


T 




torque about local x-axis 


u 




deflection in x-direction 


u 




unit dyadic 


°x 


UX 


unit vector in local x-direction 


V 


UXP 


unit vector in local x'-direction 


u „ 

x" 


UXDP 


unit vector in local x"-direction 


U 


U, U1 


transfer matrix 


V 




deflection in local y-direction 


V 




field matrix 


V 

y 




shear force in local y-direction 


V 

z 




shear force in local z-direction 
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TABLE A.l (Continued) 



TEST 


PROGRAM 




SYMBOL 


VARIABLE 


DESCRIPTION 


w 




deflection in local z-direction 


z 




state vector 


a 


AL 


x-axis coordinate transformation angle 


a l 


AL1 


x-axis coordinate transformation angle 


Y 


GAM 


z-axis coordinate transformation angle 


A 




frequency determinant 


0 




angular deflection about z-axis 


y 


UU 


mass per unit length of pipe 


6 


PR 


Poisson's ratio 


£ 


SDF 


shear distribution factor 


p 


RHOM 


mass density 






angular deflection about x-axis 






angular deflection about y-axis 


10 


W1,W2 


circular frequency 




AA 


array of section properties 




ARCL 


arc length of curved pipe 




BL 


distance from bend point to tangent 
point 




CTX12 


12 x 12 x-axis coordinate transformation 
matrix 




CTZ12 


12 x 12 z-axis coordinate transformation 
matrix 




D 


value of frequency determinant 




D1 


outside pipe diameter 




D2 


inside pipe diameter 
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TABLE A.l (Continued) 



TEXT PROGRAM 

SYMBOL VARIABLE 


DESCRIPTION 


DD 


frequency determinant array 


FN 


array of natural frequencies 


GO 


gravitational constant 


ISEC 


section identifier code 


JBC 


branch boundary condition code 


K 


boundary condition sequence number 


KB 


branch repeat discriminant 


LBC 


left boundary condition code-main member 


LBR 


branch point identifier discriminant 


M 


dimension of array of frequency determinant 
values 


MBC 


right boundary condition code-main member 


MCS 


curved subsection override discriminant 


MSR 


shear deflection/rotary inertia discriminant 


MST 


straight section test discriminant 


N 


section number 


NB 


vector of boundary condition codes 


NIB 


initial unit vector discriminant for a 
branch 


NID 


point identifier code 


NIT 


number of iterations required for Mueller's 
method 


NIV 


initial unit vector discriminant for main 
member 


NK 


dual branch discriminant 


NMO 


number of modes requested in analysis 
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TABLE A.l (Continued) 



TEXT PROGRAM 
SYMBOL VARIABLE 


DESCRIPTION 


NN 


iteration section number 


NPROB 


number of systems to be analyzed by 
program 


NSEC 


section number 


NSS 


number of curved pipe subsections 


PX 


X-coordinate of a working point 


PY 


Y-coordinate of a working point 


PZ 


Z-coordinate of a working point 


RHO 


weight density 


UA 


point matrix for a curved pipe 


UB , UCY 


field matrix for a curved pipe 


UST 


straight pipe transfer matrix 


WINCR 


frequency increment 


IVN 


natural frequency 


X 


array of frequencies for output plot 


Y 


array of frequency determinant values 
for output plot 
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STRAIGHT PIPE TRANSFER 



APPENDIX B 



CATALOG OF TRANSFER MATRICES USED IN THE PROGRAM 
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CURVED PIPE FIELD MATRIX 
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APPENDIX C 



DESCRIPTION OF PROGRAM 
1. General Remarks 

Program VIBREL is a double precision FORTRAN language 
digital computer program which is capable of performing a 
frequency analysis of a three-dimensional piping system. 

The program accomplishes this analysis by means of the 
transfer matrix method in which each section of pipe is charac- 
terized by a matrix known as a transfer matrix. These matrices 
which are functions of frequency are multiplied together 
in succession to form a transfer matrix which is characteris- 
tic of the entire system. To combine satisfactory accuracy 
with minimum computer time, a lumped mass model is used for 
curved sections of pipe while for straight pipe, a distributed 
mass approach is employed. Shear deflection and rotary 
inertia are optional in the model. 

The solution is obtained by incrementing the frequency 
and forming the system transfer matrix and frequency determi- 
nant with each increment. The sign of the frequency determi- 
nant is used to control the search for and convergence to 
system natural frequencies. The convergence acceptability 
criterion may be set arbitrarily by the user depending 
on the accuracy desired. (See section 3.b of this appendix.) 

As one would expect, the greater the accuracy, the more 
computer time required. 



The system is defined by working points which are 
read into the program with the intensive and extensive 
properties of the section of pipe preceding it. The details 
for this procedure are contained in Appendix E. The inten- 
sive and extensive properties may not vary along a single 
section of pipe but may vary from section to section. 

Projective complexities, including branches emanating from 
branches, cannot be handled directly by the program. 

Three boundary conditions, fixed, free, or ball joint, 
may be specified by coded input for the ends of the main 
member and the remote ends of the branches. Alternative 
end conditions are discussed in section 3.b of this appendix. 

2. Program Structure 
a. Main Program 

A simplified flow chart is shown in Fig. C.l. This 
discussion will pertain to that diagram. 

Data which contains working points, intensive and 
extensive properties, and boundary condition codes is read 
into the program as it is required for computation. Local - 
coordinate systems are set up at the starting point, prior 
to and following each directional change, and at branch 
junction points. From the working point locations and unit 
vectors describing the local coordinate systems, section 
lengths and coordinate rotation angles are computed. 

For computational purposes piping bends are sub- 
sectioned based on the length-to-diameter ratio in the following 
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Figure C.l: Simplified Flow Diagram for Main Program 
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manner : 



0 < L/D < 1 1 SUBSECTION 

1 < L/D <_ 3 2 SUBSECTIONS 

3 < L/D < 6 3 SUBSECTIONS 

If L/D is greater than six, the number of subsections is 

computed as twice the number of mode frequencies sought, up 
to a maximum of twelve. The user may, however, choose to 
override the number of subsections as determined by the 
program through the use of the parameter MCS described in 
Appendix E. 

As the lengths, subsections, and rotation angles 
are computed, they are stored, together with intensive and 
extensive properties, in an array. This array is passed 
via common storage (dotted line in Fig. C.l) to function 
subroutine FRDET (FRequency DETerminant) which later forms 
and evaluates the frequency determinant for a given value 
of frequency. 

When the geometry for the entire system has been 
computed, frequency iteration begins. The starting frequency 
is arbitrarily set at 0.1 radians per second. To compute 
the starting frequency increment, the program constructs 
a synthetic straight pipe, with a length equal to the 
total length of the main member, and having outside diameter, 
wall thickness, and intensive properties equal in magnitude 
to a weighted average of the main member sections. One 
sixth of the fundamental frequency of this synthetic pipe 



with fixed ends is taken as the starting frequency 
increment. The system transfer matrix and frequency deter- 
minant are formed by function subroutine FRDET for each 
new frequency. When a sign change is detected in the 
frequency determinant, an iteration process using Mueller's 
method of successive bisections and inverse parabolic inter- 
polation is used for convergence to the natural frequency. 

The convergence acceptability criterion is taken as the start- 
ing frequency increment divided by ten thousand. 

After a root has been located, the search begins 
anew with the starting frequency equal to the value of the 
natural frequency plus one radian per second. When the second 
mode frequency has been located, the frequency increment is 
changed to one tenth the difference between the previous 
two natural frequencies. 

The search process continues until the specified 
number of frequencies has been found or until the significant 
figure capacity of the computer has been exceeded. Then the 
mode frequencies are printed and the pairs of frequencies 
and frequency determinant values used in the search are 
sorted and plotted, 
b. Subroutines 

The program utilizes sixteen subroutines in the 
process of computation. A brief description follows: 

ANGLE Computes the sign and the angular 

separation of two sets of unit vectors 
which have a common axis. 
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COORDX 


Forms an n x n (n is some multiple of 
3 and < 12) coordinate transformation 
matrix for a rotation about a local 
x-axis . 


COORDZ 


Forms an n x n coordinate transformation 
matrix for a rotation about a local z 
axis . 


CURMAT 


Constructs the field matrix for a curved 
section of pipe. 


DETER 


Computes the value of an n x n determi- 
nant; used to evaluate the frequency 
determinant . 


DROOT 


Performs successive bisections and inverse 
parabolic interpolation to locate the 
zeros of the frequency determinant after 
a sign change has been detected. 


INVERT 


Inverts the matrix R^ to be used in 
computing the branch point matrix. 


MATMUL 


Multiplies two matrices together. One 
of the matrices must be square. 


POINT 


Constructs a point matrix for a curved 
section of pipe. 


PRPLOT 


Takes the array of frequencies and 
frequency determinant values calculated 
during iteration and plots them length- 
wise with the printer. Scaling is done 
automatically and the length of plot 
may be changed to suit the user. 


SORT 


Sorts the array of frequencies and 
frequency determinant values from 
iteration for plotting on the printer. 


STATEM 


Constructs the state matrix for the 
appropriate fixed, free, or ball joint 
boundary conditions. 


STRMAT 


Constructs the transfer matrix for a 
straight section of pipe. 


UVEC 


Computes a set of orthogonal unit vectors 
from two pipe section vectors which 
represent a directional change in the 
piping . 
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IVSTART Computes the value of the frequency 

determinant for a previously constructed 
synthetic pipe equal in length to the main 
member and having weighted average values 
of the properties of the main member 
sections. This frequency value is used 
in selecting the starting frequency 
increment . 

FRDET For a particular value of frequency, 

this subroutine assembles the state 
matrices, transfer matrices for each 
section of pipe, and the coordinate 
transformation matrices. It then 
multiplies them together in sequence 
and calculates the frequency determinant. 

An interesting phenomenon occurs in the formation of 
the frequency determinant for a single straight pipe. 

Because of the symmetry of a pipe in the y and z directions, 
the two flexural submatrices in the straight pipe transfer 
matrix are identical, which causes the frequency determinant 
to maintain the same sign on either side of the flexural 
natural frequencies. This isolated case is handled in sub- 
routine STRMAT through the use of the straight section test 
discriminant and appropriate exclusion of one of the flexural 
submatrices . 



3. Remarks on Instructions for Program Use 
a. General Remarks 

Detailed instructions for program use are contained 
in the listing of Appendix E. It was intended for these 
instructions to be part of the program in order that any 
user, given only the program listing or the computer card 
deck, could employ the VIBREL package with no prior knowledge 
of transfer matrices and without perusal of this thesis. 



63 



Some clarification of the instructions for program 
use is justified for cases where properties or geometry 
change within a piping bend. First, consider the simple 
section of piping shown in Fig. C.2. Assuming that no change 





o 



Figure C.2: Piping Bend with No Change in Properties 



in cross section or properties occurs between points a and 
d, the only points which must be listed as working points 
for program input are a, o, and d. Using these points and 
the radius of curvature, program VIBREL will calculate 
lengths ab and cd as well as the arc length and angle of the 
bend . 

Now consider the section of piping of Fig. C.3 where 
the cross section or properties change at point c but remain 



Figure C.3: 



a 
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constant between points a and c and between points c and e. 

In this instance points a, o^, O 2 , and e must be listed as 
working points for the program. VIBREL calculates straight 
lengths ab and de as well as the information required for 
curved lengths be and cd. It should be obvious from this 
example that for the general case of two curved sections occur- 
ring sequentially that the intersection point (point c in 
Fig. C.3) need not be listed. 

If the geometry input is such that negative lengths 
are computed, an error message is printed. One half inch 
leeway is allowed for negative lengths to take into account 
any errors which may occur in measurement, 
b. Program Flexibilities 

(1) Approximating Smal-1 Piping Accessories 

Small valves, flanges, couplings and other piping 
accessories whose center of mass is relatively close to the 
centroidal axis of the pipe can be approximated using a 
straight section of pipe. The length, mass density, diameter 
and wall thickness can be artificially varied to give the 
mass of the item. The resulting variation in stiffness 
between an actual flange and artificial pipe, for example, 
usually will not have a significant effect on the results 
because of its short length. 

(2) Approximating Alternative End Conditions 

When the ideal boundary conditions are applied 
to the state vector on the end of a system, it must have 
six zero and six non-zero elements. Actual field conditions 
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may dictate that this is not the case because of flexible 
mountings or connection to machinery. A good approximation 
in these circumstances may be made by attaching one end of 
an artificial straight pipe to the system piping keeping its 
other end fixed. The mass of the artificial pipe should be 
small with the extensive properties and elastic modulus 
chosen to approximate the stiffness of the mounting. 

The program can then use the ideal boundary 
conditions with little loss of accuracy from the real situation. 

(3) Program Modifications 

In some instances, a user may desire to change 
such quantities as frequency increment or input format to 
suit his particular needs. For this reason, Table C.l 
lists the more important program variables and what card or 
cards to make these changes to. 



TABLE C.l 



Modification 


Card 

Numbers 


Variables 

Affected 


Maximum number of points for 
print plot 


2890 


x (400) , y (400) 


Maximum number of piping sections 
in analysis 


2910 


AA (100,10) 


Maximum number of branches in 
analysis 


2920 


NB ( 50 ) 


Input format 


3210-50 




Output format (Input Data) 


3390-420 




Output format (Properties and 
Geometry) 


7750-60 




Output format (Natural 


8740 





Frequencies) 
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TABLE C.l (Continued) 



Modification 


Card 

Number 


Variables 

Affected 


Starting frequency increment 


86S0 


WINCR 


Starting frequency 


8760 


W1 


Com^ergence acceptability 
criterion 


9120 


EPS 


Maximum iterations to 
convergence in DROOT 


3670 


NIT 


Starting Frequency after 
location of natural frequency 


9330 


W1 


Frequency increment after 
second mode 


9320 


WINCR 


Number of lines of print plot 


14730 


RY 
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APPENDIX D 



PROGRAM FLOW DIAGRAMS 



1. General Remarks 

Flow diagrams are included for the main program and 
function subroutine FRDET. Subroutines DETER, DROOT, 

INVERT, and PRPLOT are standard package subroutines with 
minor modifications made to accommodate them to program 
VIBREL. Since their equivalent or an alternative routine 
may be substututed for them, flow diagrams are not included. 
Adequate comment statements are, however, contained in these 
subroutines to define their structure. Other subroutines 
for which no flow diagrams appear, incorporate, at most, two 
decisions in their structure* - 

2. Main Program Flow Diagram 
Pages 6? through 76. 

3. Function Subroutine FRDET Flow Diagram 
Pages 77 through 79 . 
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APPENDIX E 



PROGRAM LISTING 



1. Index 



SECTION 


SEQUENCE NUMBERS 


PAGE 


Instructions for Use 


10-2810 


81 


Main Program 


2820-9650 


86 


Subroutines 






ANGLE 


9660-9860 


101 


COORDX 


9870-10250 


101 


COORDZ 


10260-10640 


102 


CURMAT 


10650-11700 


103 


DETER 


11710-12060 


105 


DROOT 


12070-13200 


106 


INVERT 


13210-13980 


108 


MATMUL 


13990-14130 


110 


POINT 


14140-14570 


110 


PRPLOT 


14580-15910 


111 


SORT 


15920-16150 


114 


STATEM 


16160-16370 


115 


STRMAT 


16380-17370 


115 


UVEC 


17380-17630 


117 


WSTART 


17640-17830 


118 


FRDET 


17840-20090 


118 
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APPENDIX F 

PROGRAM ACCURACY AND INTEGRITY ANALYSIS 
1. General Remarks 

This appendix tabulates by mode frequency (radians/second) 
the results of the accuracy and integrity analyses of various 
piping systems from simple to complex. 

The comparison frequency values for the simple straight 
section systems were obtained by solution of the differential 
equations for longitudinal, torsional and flexural vibrations 
with application of the appropriate boundary conditions. These 
solutions were evaluated on the digital computer with accuracy 
to at least as many places as listed in the tabulation. Evalu- 
ation of the comparison solution and VIBREL solution were 

2 

carried out using a value of 32.174 FT/SEC for the gravita- 
tional constant and a value of tt to 15 significant figures. 

Comparison values for the curved pipe analysis were obtained 
from three sources: Fink [5], Kim [6], and Ojalvo [8], Those 

extracted directly from [5] and [6] were compiled using a lumped 
mass model and digital computation. From graphs and formulas 
(based on the theory of clamped ring segments) appearing in [8], 
the other natural frequency values for the curved pipe were 
calculated . 

Several comparison frequencies for a single branch system 
were obtained from [6] . 

Shear deflection and rotary inertia were neglected except 
where noted. 
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2. SYSTEMS ANALYSIS 



a. STRAIGHT PIPE WITH VARIOUS END CONDITIONS 







END CONDITIONS: 
FIXED-FIXED 



SYSTEM PARAMETERS: 



DIAMETER (IN) 2.0 

WALL THICKNESS (IN) .125 

LENGTH (IN) 200. 

MODULUS OF ELASTICITY (PSI) 30 x 10 

POISSON'S RATIO - .25 

WEIGHT DENSITY (LB/FT ) 490. 



COMPARISON RESULTS: 



MODE 


ANALYTICAL 


VI BREL 


1 DIFFERENCE 


1ST 


FLEXURAL 


75.10461420 


75.10461590 


.00000 


2ND 


FLEXURAL 


207.02876110 


207.02876110 


. 00000 


3RD 


FLEXURAL 


405 .85934940 


405.85914210 


.00000 


4TH 


FLEXURAL 


670.90578980 


670 . 90577600 


. 00000 


5TH 


FLEXURAL 


1002.21750380 


1002 . 21750410 


.00000 


6TH 


FLEXURAL 


1399.79137910 


1399.79140510 


. 00000 


7TH 


FLEXURAL 


1863.62757630 


1863.62799710 


.00002 


1ST 


TORSIONAL 


2007.8325054 


2007.8275305 


. 00025 


8TH 


FLEXURAL 


2393.7260868 


2393.7157896 


.00043 


9TII 


FLEXURAL 


2990.0869112 


2990.3653337 


.00931 


1ST 


LONGITUDINAL 


3174.6619386 


3174.6560265 


.00019 



DUAL ANALYSIS 


RESULTS: 








STARTING AT 


STARTING AT 




MODE 


LEFT END 


RIGHT END 


% DIFFERENCE 


1 


75.1046142 


75 . 1046142 


. 00000 


2 


207.0287611 


207.0287611 


.00000 


3 


405.8591494 


405 .8591494 


.00000 


4 


670.9057898 


670.9057898 


. 00000 


5 


1002.2175038 


1002 .2175038 


.00000 


6 


1399 . 7913791 


1399.7913791 


.00000 


7 


1863.6275763 


1863.6275763 


.00000 
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END CONDITIONS: 
FIXED-FREE 



SYSTEM PARAMETERS: 
SAME AS (1) 




COMPARISON RESULTS: 



MODE 


ANALYTICAL 


VIRREL 


% DIFFERENCE 


1 


11.8028657 


11.8028696 


.00003 


2 


73.9673211 


73.9673092 


.00002 


3 


207.1106480 


207.1106397 


. 00000 


4 


405.8541957 


405 .8541953 


.00000 


5 


670.9060651 


670.9060651 


. 00000 


6 


1399.7913799 


1399 . 7913885 


.00000 



DUAL ANALYSIS 


RESULTS: 








STARTING AT 


STARTING AT 




MODE 


LEFT END 


RIGHT END 


% DIFFERENCE 


1 


11.8028696 


11 .8028696 


.00000 


2 


73.9673092 


73.9673092 


.00000 


3 


207 .1106397 


207.1106397 


. 00000 


4 


405.8541953 


405 .8541953 


.00000 


5 


670.9060651 


670.9060651 


. 00000 


6 


1399.7913885 


1399.7913885 


.00000 
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SYSTEM PARAMETERS: 




END CONDITIONS: 
FIXED-PINNED 



COMPARISON RESULTS: 



MODE 


ANALYTICAL 


VIBREL 


% DIFFERENCE 


1 


51.7571905 


51.7571773 


. 00003 


2 


167 .7264474 


167 .7264180 


.00002 


3 


349.9478448 


349.9478460 


.00000 


4 


598.4315201 


598.4315226 


.00000 


5 


913.1775124 


913.17753170 


.00000 


6 


1003.9162520 


1003.9164796 


.00002 



DUAL ANALYSIS 


RESULTS: 








STARTING AT 


STARTING AT 




MODE 


LEFT END 


RIGHT END 


% DIFFERENCE 


1 


51 .7571773 


51.7571773 


.00000 


2 


167.7264180 


167.7264180 


.00000 


3 


349.9478460 


349 . 9478460 


.00000 


4 


598.4315226 


598.4315226 


. 00000 


5 


913.17753170 


913.17753170 


.00000 


6 


1003.9164796 


1003.9164796 


.00000 
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Vu ViU V Lil/ 



1 X J 1 , 




END CONDITIONS: FIXED- FIXED 

SYSTEM PARAMETERS: 



PIPE DIAMETER (IN) 2.0 

WALL THICKNESS (IN) .125 

INCLUDED ANGLE OF ARC (DEG) 270 

MODULUS OF ELASTICITY (PSI) 30 x 10 6 

POISSON'S RATIO .25 

WEIGHT DENSITY (LB/FT -5 ) 556 



COMPARISON RESULTS 



a . MODE 1 








SOURCE 


FREQUENCY 


% 


VARIATION WITH 
VIBREL VALUE 


VIBREL 


37.2234058 






OJALVO [8] 


38.5930983 




3.68 


FINK [5] 


38.8773520 




4.44 


KIM [6] 


35.3639740 




-5.00 


b . MODE 2 








SOURCE 


FREQUENCY 


1 


VARIATION WITH 
VIBREL VALUE 


VIBREL 


72.1914615 






OJALVO 


69.9962862 




-3.04 


FINK 


73.3297307 




1 . 58 


KIM 


66.6993490 




-7.61 
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c. BRANCHED SYSTEMS 

(1) SINGLE BRANCH SYSTEM 




SYSTEM PARAMETERS (ALL 3 SFCTIONS): 



DIAMETER (IN) 2.0 

WALL THICKNESS (IN) .125 

LENGTH (IN) 100. 

MODULUS OF ELASTICITY (PSI) 30 x 10 6 

POISSON'S RATIO .25 

WEIGHT DENSITY (LB/FT 3 ) 556. 



COMPARISON RESULTS: 



MODE 


KIM [6] 


VI BREL 


% DIFFERENCE 


2 


194.25470558 


194 . 2888499 


.01758 


5 


278.61078241 


278.6390155 


. 01013 


8 


628 . 59556355 


628.9417782 


.05508 



DUAL ANALYSIS 


RESULTS: 








STARTING AT 


STARTING AT 




MODE 


LEFT END 


RIGHT END 


% DIFFERENCE 


1 


67.6456627 


67.6456627 


. 00000 


2 


194.2888474 


194 . 2888474 


.00000 


3 


198.0817333 


198 . 0817333 


. 00000 


4 


242 . 0773169 


242 . 0773169 


. 00000 


5 


278.6390155 


278.6390155 


.00000 


6 


281.9864461 


281.9864461 


. 00000 


7 


370.7949928 


370 . 7949928 


.00000 


8 


628.9417782 


628.9417782 


.00000 


9 


634.7478785 


634.7478786 


. 00000 


10 


689.7373111 


689.73731 11 


. 00000 
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(2) DUAL BRANCH SYSTEM 




END CONDITIONS: 

MAIN MEMBER 
BRANCH #1 
BRANCH #2 



- FIXED-FIXED 

- PINNED 

- FREE 



SYSTEM PARAMETERS: 



PARAMETER 



SECTION 

1 



SECTION SECTION 

2 3 



SECTION 

4 



DIAMETER (IN) 

WALL THICKNESS (IN) 
LENGTH (IN) 

MOD. OF ELAST. (PSI) 
POISSON'S RATIO 
WT. DENSITY (LB/FT 3 ) 



2.0 

.125 

80.0 
30 x 10 
.30 

490.0 



2.0 

. 12 5 

50.0 

30 x 10 
.30 

490.0 



2.0 

.125 

31.6 

30 x 10 

.30 

490.0 



2.0 

.125 

60.0 

30x10 

.30 

490.0 



DUAL ANALYSIS 



MODE 

1 

2 

3 

4 

5 

6 



SULTS* : 

STARTING AT 
LEFT END 

110.5015347 

276.6928874 

314.3597336 

417.0897739 

463.5288123 

675.6526889 



STARTING AT 
RIGHT END 

110.5015347 

276.6928874 

314.3597336 

417.0897739 

463.5288123 

675.6526889 



% DIFFERENCE 

.00000 
.00000 
.00000 
.00000 
. 00000 
.00000 



* SHEAR DEFLECTION AND ROTARY INERTIA ARE CONSIDERED IN THIS 
ANALYSIS. 

+ COORDINATES SHOWN IN DIAGRAM ARE IN INCHES. 
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d. 



COMPLEX SYSTEM 

END CONDITIONS: 



SYSTEM PARAMETERS 
(ALL SECTIONS) : 



MAIN MEMBER - FIXED- FIXED 
BRANCH - FIXED 

Os Q,ZS,-ZS') 



DIAMETER (IN) 

WALL THICKNESS (IN) 



2.25 

.035 



ELASTIC MODULUS (PSI) 30 x 10 

. 30 
491.5 




DUAL ANALYSIS 



• MODE 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 



RESULTS: 

STARTING AT 
LEFT END 

73.0350449 

102.7404907 

142.5065200 

201.7233189 

239.1786608 

327.8210841 

331.8165320 

364.5803428 

428.5622381 

520.1734071 

542.4026997 

565.0762600 

574.0701684 

648.8240390 

754.8640271 

1052.0436077 

1210 .9202599 

1250 .1288524 



STARTING AT 
RIGHT END 

73.7545503 

102.9266533 

146.3168530 

201.3564078 

239.1641100 

321.1763486 

337.3328726 

364.1833110 

429.1841008 

521.2514969 

544.3290021 

564.3679439 

571.0709491 

651.2831922 

759.4399716 

1052.3129250 

1210.7321375 

1247.3298693 



% DIFFERENCE 

.985 
.181 
2.674 
.182 
.006 
2.027 
1.662 
.109 
. 145 
.207 
.355 
.125 
.522 
.379 
.606 
.026 
.016 
.224 
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APPENDIX G 



SAMPLE PROBLEM 



1. Piping System 

The piping used in this example is identical to the 
complex system shown in Appendix F. It includes several corners 
and bends and has one branch junction point. The boundary 
conditions are fixed for the two ends of the main member and 
the remote branch end. 

2. Data Cards 

Shown in Fig. G.l is the arrangement of data cards 
as they would appear following the program deck. Note that 
Card #1 specifies one system to be analyzed, and Card #2 
asks for four mode frequencies in the analysis. Since the 
curved subsection override is not utilized, the program will 
go ahead and determine the number of curved pipe subdivisions 
according to the length- to- diameter ratio. Shear deflection 
and rotary inertia are not included in this problem analysis. 
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Figure G.l: Data Deck for Sample Program 
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************************************************************************* 



♦ PROGRAM VIBREL 

* 

*************************************** 



C.O.RUOOLF * 

* 

********************************** 



PROGRAM VIBREL COMPUTES THE NATURAL FREQUENCIES OF VI8RATION 
OF ANY RANOOMLY ARRANGEO 3-0 1 M t NS IONAL PIPING SYSTEM 
PIPING HANGERS ANO PROJECTIVE COMPLEXITIES ARE NOT INCLUOEO 



***************** ********************************************* ************ 

PROBLEM n 1 

************************************************************************** 

THE EFFECTS OF SHEAR DEFLECTION ANO ROTARY INERTIA ARE NOT CONSIOEREO IN THIS PROBLEM 
******* ******************************************************************* 

INPUT OATA: 





X 


Y 


Z 


RAD OF 


ELASTIC 


POISSONS 




WALL 




LEFT BC 


RIGHT BC 


10 


COORO 


COORO 


COCRO 


CUR V 


MOOULUS 


RATIO 


o 

o 


THICK 


WT/VOL 


COOE 


COOE 


1 


0.0 


0 . 0 


0.0 














1 


1 


3 


20.0 


0.0 


0.0 


0.0 


0.3000 OB 


.30 


2.25 


0.0350 


491.5 






5 


35.0 


0.0 


0.0 


0.0 


0.0 


.0 


0.0 


0.0 


0.0 






5 


60.0 


25.0 


-25.0 


0.0 


0.0 


.0 


0.0 


0.0 


0.0 






3 


60.0 


25.0 


10.0 


0.0 


0.0 


.0 


0.0 


0.0 


0.0 






A 


60.0 


25.0 


40.0 


3.0 


0.0 


.0 


0.0 


0.0 


0.0 






6 


90.0 


-9.0 


55.0 


0.0 


0 .0 


.0 


0.0 


0.0 


0.0 






2 


140.0 


20.0 


-10.0 














1 




4 


140. 0 


20.0 


20.0 


3.0 


0.0 


.0 


0.0 


0.0 


0.0 






5 


1 10. 0 


20.0 


20.0 


0.0 


0.0 


.0 


0.0 


0.0 


0.0 






B 


90.0 


-9.0 


55.0 


0.0 


0.0 


.0 


0.0 


0.0 


0.0 






B 


100.0 


-20.0 


60.0 


0.0 


0.0 


.0 


0.0 


0.0 


0.0 







************************************************************************** 



PROPERTIES ANO GEOMETRY: 



SECT 

NR 


I 0 


ELAST IC 
MOOULUS 


POISSONS 
RAT 10 


WT/VOL 


0 . 0 . 


WALL 

THICK 


LENGTH OR 
RAO OF CUR V 


ALPHA1 


GAMMA 


NR OF SUBSEC 
OR ALPHA 


1 


1.0 


0.3000 


OB 


.30 


491.5 


2.25 


0.0350 


20.000 


0.0 


0.0 


0.0 


2 


3 .0 


0.3000 


08 


.30 


491 . 5 


2.25 


0.0350 


15.000 


-2.356 


0.955 


0.0 


3 


3.0 


0.3000 


08 


.30 


491.5 


2.25 


0.0350 


43.301 


2.094 


2. 1B6 


0.0 


4 


1.0 


0.300D 


08 


.30 


491.5 


2.25 


0.0350 


35.000 


0.0 


0.0 


0.0 


5 


1 .0 


0.3000 


OB 


.30 


491 .5 


2.25 


0.0350 


27. B33 


0.0 


0.0 


0.0 


6 


2.0 


0.3000 


06 


.30 


491.5 


2.25 


0.0350 


3.000 


1.50B 


1.251 


2.000 


7 


5.0 


0.300D 


08 


.30 


491 . 5 


2.25 


0. 0350 


45.592 


0.0 


0.0 


0.0 


B 


1.0 


0.3000 


08 


.30 


491 . 5 


2.25 


0.0350 


27.000 


0.0 


0.0 


0.0 


9 


2.0 


0.300D 


OB 


.30 


491.5 


2.25 


0.0350 


3.000 


0.0 


1.571 


2.000 


10 


3.0 


0.3000 


OB 


.30 


491 .5 


2.25 


0.0350 


27.000 


2.450 


1 .156 


0.0 


1 1 


7.0 


0.?00D 


08 


.30 


491.5 


2.25 


0.0350 


49.659 


-0. 372 


-1.177 


-0. B47 


12 


4.0 


0.3000 


08 


.30 


491.5 


2.25 


0.0350 


15 .6B4 


0.0 


0.0 


0.0 



************************************************************************** 
BOUNOARY CONOITION COCES: 



LEFT B.C. COOE = 1 
RIGHT B.C. COOE = 1 
BRANCH # 1 B.C. COOE= 1 



************************************************************************** 



PIPING SYSTEM NATURAL FREQUENCIES: 



MOOE 


FREQUENCYIRAO/SEC) 


FREQUENCY DETERMINANT 


1 


73.0350449 


0. 122780-22 


2 


102.7404907 


-. 369030-21 


3 


142.5065200 


0. B5913D-21 


4 


201.7233189 


0.524240-22 



t 
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GRAPH OF THE value OF FREQUENCY OETERMINANT VS. FREQUENCY 



-5.113E-1? -2.538E-17 3.585E-19 2.610E-17 5.184E-17 

0.0 ♦ • ♦ 



1.270E 01 



2.539E 01 



3.809E 01 



5.079E 01 



6.3*9E 01 



7.618E 01 



8.888E 01 



1.01 6E 02 



1.143E 02 



1.270E 02 



1.397E 02 



1.52 4E 02 



1 • 65 IE 02 



1.778E 02 



1.905E 02 



2.032E 02 ♦ I ♦ 

-5.113E-17 -2.538E-17 * 3.585E-19* * 2.610E-17 * 5.18AE-17 



0.0 



1.270E 01 



2.539E 01 



3.809E 01 



5.079E 01 



6.3A9E 01 



7.618E 01 



1.016E 02 



1.1A3E 02 



1.270E 02 



1.397E 02 



1.52AE 02 



1.651E 02 



1.778E 02 



1.9 05E 02 



2 . 032 E 02 
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